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Using molecular-dynamics simulation we simulate nanoindentation into the three principal sur- 
faces - the (100), (110) and (111) surface - of Cu and Al. In the elastic regime, the simulation 
data agree fairly well with the linear elastic theory of indentation into an elastically anisotropic 
substrate. With increasing indentation, the effect of pressure hardening becomes visible. When 
the critical stress for dislocation nucleation is reached, even the elastically isotropic Al shows a 
strong dependence of the surface orientation on the force-displacement curves. After the load drop, 
when plasticity has set in, the influence of the surface orientation is lost, and the contact pressure 
(hardness) becomes independent of the surface orientation. 
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I. INTRODUCTION 

In a seminal paper, Hertz investigated the elastic interaction between a sphere of radius R with an elastic isotropic 
solid. 1 He found that the force F acting on the sphere normally to the surface and the displacement d into the surface 
are related by the so-called Hertz law 

F = ^Md 3/2 y/R. (1) 

In this relation, a single materials parameter, the so-called indentation modulus M describes the materials clastic 
response. 1-4 For a rigid indenter, it may be expressed in terms of the Young's modulus E and the Poisson ratio v of 
the substrate as 

Hertz also determined the contact pressure; it is defined by the ratio of the normal force F divided by the contact 
area projected into the surface plane. Hertz obtained 




P=^Z M \-5- (3) 



We note that Hertz' analysis is exact in the limit of d <C R. 

Nowadays indentation experiments can be performed into single crystals. 5 ' 6 The question then arises in how far 
Hertz' analysis also describes the indentation of crystalline - and thus by definition anisotropic - materials. An 
analytical extension of Hertz' analysis to anisotropic materials is, however, non-trivial. Willis ,s appears to have gone 
farthest in the analysis; however, no analytical results are available, and one has to resort to numerical procedures even 
in the simplest case, i.e., for transversally isotropic media. Vlassak and Nix 9 ' 10 evaluated these results numerically 
for the specific problems of a flat circular punch and an axisymmetric paraboloid indenter; the latter approximates a 
spherical indenter for small indentation depths. They showed that in this case the Hertzian law Eq. (1) holds with 
a modified indentation modulus M(Q) which depends on the surface orientation $1, and calculated this quantity 
numerically. In their experiments on Au single crystals, Kiely and Houston found considerable deviations between 
the experimental results on a Au substrate and the theoretical predictions. 5 Recently, Tsuru and Shibutani employed 
molecular-dynamics (MD) simulations to study indentation into an anisotropic half-space, and demonstrated that 
indentation into fee crystals indeed depends on the surface orientation. 11 ' 12 Other authors also investigated the 



2 



question in how far the homogeneous nuclcation of dislocations under the indenter, i.e., the onset of plasticity, is 
influenced by the surface orientation. 13-15 

In the present paper, we study these questions for two fee materials, Al and Cu, which show a widely differing 
anisotropy. This allows us to discuss quantitatively in how far the Hertzian law, (1), is obeyed. Finally, we demonstrate 
that with increasing penetration into the target, the indenter measures the bulk hardness of the system and the surface 
orientation loses its influence. 



II. METHOD 



A. Simulation 



We chose two different fee materials, Al and Cu, for our study, which differ strongly in their elastic anisotropy. 
Many-body potentials of the embedded-atom type are used to model these metals; both the Al potential 16 and the Cu 
potential 17 reproduce the zero-temperature clastic moduli given in Table I. The fee substrate has approximately cubic 
shape with side lengths of around 25 nm; it contains roughly 1.35 x 10 6 atoms. We checked in a series of simulations 
that our crystallite size is large enough to obtain reliable results for the indentation process. We found that a careful 
relaxation of the crystal before starting the indentation process to Pij < 10 -5 GPa and temperatures « 1 K was 
mandatory to obtain reliable and reproducible results. Lateral periodic boundary conditions have been applied. At 
the bottom, atoms in a layer of the width r cut have been constrained to i^normai = 0. 

The indenter is modelled as a repulsive sphere. We chose a non-atomistic representation of the indenter. Its 
interaction potential with the substrate atoms is described by 18 



= (4) 

The indenter radius was set to R = 8 nm, and the indenter stiffness to k ~ 3 eV/A 3 . We checked that our results 
are only weakly influenced by the exact value of the contact stiffness, as long as it is in the range of 1 - 10 eV/A 3 . 

The simulations have been performed using a modified version of the LAMMPS code, 19 using the so-called 
displacement-controlled approach. 20,21 The indenter is advanced every At = 2 ps by a fixed amount of 5 = 0.256 
A (<C lattice constant) instantaneously, corresponding to an average indentation speed of v = 12.8 m/s. The sub- 
strate then relaxes for the ensuing time of At to the new indenter position. 



B. Elastic properties 



The materials considered here crystallize in the cubic fee structure. Their elastic behaviour is therefore completely 
described by their elastic constants, cu, C12, and C44. Their zero-temperature values, are given in Table I. These two 
materials differ in particular in their elastic anisotropy, defined as 



X = (5) 

en - C12 

which measures both the orientation dependence of the elastic modulus and of the shear moduli. As Table I shows, 
Cu is quite anisotropic, X = 3.22, while Al is nearly isotropic, X = 1.20. Fig. 1 demonstrates the anisotropy of the 
two materials by plotting Young's modulus, E(fl), in its dependence on orientation f2. It is given by 



where E is the orientation-averaged Young's modulus, which is related to the average shear modulus G (cf. Eq. (10)) 
and the average Poisson ratio v (cf. Eq. (11)) via 



E = 2G(v + l), 



(7) 



and 
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r(n) = Oi0 2 + ^2 JV + o 3 (8) 

describes the angular dependence. 22 

For reference purposes, it is interesting to define pertinent isotropic properties, such as they would apply to a 
polycrystalline sample: 9 the bulk modulus B, 

B=^±^l, (9) 
an average shear modulus G, defined as the arithmetic mean over the Voigt- and Reuss-averaged shear moduli 

q_\( Cll - C12 + 3c 44 5(Cn - Ci 2 )c 4 4 \ 

~2^ 2(G + 3B) 4c 4 4 + 3(cii -da)/ ' 
and an average Poisson ratio v as appropriate for an isotropic solid 

3B-2G 

Indentation moduli M(f!) for the three crystal orientations have been determined numerically by Vlassak and 
Nix. 9 ' 10 We evaluate these moduli for Al and Cu and give the results, denoted as M theo in Table II. The orientation 
dependence amounts to only 1 % for Al, while for Cu, deviations of up to 10 % from the isotropic value, M lso , show 
up. In agreement with the orientation dependence of the Young's modulus, Fig. 1, the (111) plane is the stiffest, and 
the (100) plane the softest; the (110) plane is intermediate, and its modulus is quite close to the isotropic value. 

Other authors 12 attempt a different procedure and define 



Thus they use the definition of the isotropic indentation modulus Eq. (2), and replace Young's modulus E by the 
orientation dependent Young modulus E(Sl). Analogously, v is replaced by the Poisson ratio f(f2), which characterizes 
the contraction transverse to ft; for the (100) and (111) orientation, this value is uniquely defined, while for the (110) 
direction a Voigt average is applied. The resulting values of this approximated indentation modulus are displayed in 
Table II. They show a considerably larger anisotropy than the Vlassak-Nix moduli; the deviation from the isotropic 
modulus vary between —9 and +6 % for Al and between —44 and +45 % for Cu. Evidently, M appr is not a good 
approximation to the true theoretical value, M thco . 



III. RESULTS 
A. Elastic regime 



In Fig. 2 we show the force-displacement curves for Al and Cu indented on the three principal surface planes. In the 
elastic regime, which is characterized by the d 3 ! 2 dependence of the force and terminates with a more or less drastic 
force maximum at a depth of around 6 — 8 A, Cu is stiffcr than Al, corresponding to the larger elastic moduli of Cu. 
Both materials show a dependence on the surface orientation, such that the stiffness is highest for the (111) plane 
and smallest for the (100) plane, in agreement with the orientation dependence of the indentation modulus, cf. Table 
II. The orientation dependence is most pronounced for Cu, corresponding to the larger anisotropy X of this material. 

We fit the elastic part of the force-indentation curve to the Hertzian law, Eq. (1), using M as a fit parameter. We 
perform two fits: (i) over the full elastic part of the F(d) curve up to d = 5 A, giving a fit parameter M MD,iul1 ; (ii) 
a fit only the initial part, d < 2 A, giving 7^f MD > mlt . The fit values are included in table II. We see that the full fit 
gives consistently larger fitted indentation moduli than the initial fit, jVf MD > fuU > jyMD.imt^ This is a sign of the onset 
of nonlinear elasticity: under pressure, a material shows an increased stiffness; this feature has been termed pressure 
hardening. 23 ' 24 This effect of pressure hardening has been demonstrated previously by comparing FEM simulations of 
the indentation process using linear vs nonlinear elasticity. 23 The fact that the indenter is not exactly a hard sphere, 
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but for numerical reasons - has to be defined with a finite stiffness in form of the polynomial indenter potential, 
Eq. (4), also causes the initial stiffness to be systematically smaller than the stiffness at larger penetration depths; 
this feature enhances the effect of the pressure hardening. 

A comparison of the simulation results with theory hence must be based on the initial fit values M MD > lnlt . Table II 
demonstrates a fairly good agreement between theory and simulation; the largest deviations amount to 14 %. We note 
that - with the exception of the Cu (100) surface - all fit values are above the theoretical values; this points at the 
possibility that even for the small fit regime of < d < 2 A, already some departures from linear elasticity may show 
up. Among other reasons which may contribute to the deviations we mention the atomistic nature of the indentation 
process, which for the indenter radius of 8 nm is not fully captured by continuum elasticity, and also the numerical 
problem of fitting the molecular-dynamics data to the Hertzian law - note that besides the elastic deformation also 
a finite offset in the displacement has to be fitted. 

In Fig. 3 we give a graphical representation of the comparison between MD simulation results, the fit curves, and the 
theoretical prediction for the case of the (100) surfaces. Both for Al and Cu, the MD data are almost indistinguishable 
from the full fit. It is also evident that deviations between the full fit and the initial fit become sizable only for d > 2 
A, as it must be. The approximate indentation moduli, Eq. (12), severely underestimate the simulation data, while 
the theoretical prediction by Vlassak and Nix gives a fair representation of the simulation data. 

B. Plasticity 

The onset of plasticity is visible in the force-displacement curves of Fig. 2 by the pronounced load drop appearing 
for the (111) and - to a lesser extent - for the (110) surface. In the following we denote the stress at which this 
load drop appears as the critical stress and the pertinent indentation as the critical indentation depth. For the (100) 
surface, dislocation nucleation sets in considerably earlier and in a more continuous way, such that any load drop is 
blurred out. 

These principal differences between the (100) and (111) surfaces and the abrupt onset of plasticity for the (111) 
surface can be explained as follows: (i) the primary glide systems (110) {111} are located at quite oblique angles 
to the di rectio n of the indentation force acting normally to the surface; the corresponding Schmid factor is only 
sm = \/2/27 = 0.27. For the (100) surface these glide systems are more easily activated, since s W o = 1.5sm = 
1/V6 = 0.41. (ii) When finally the critical indentation depth has been reached, a considerable elastic energy has built 
up due to the high stiffness of this surface. Then, upon dislocation nucleation, a stronger dislocation avalanche and 
consequently a larger plastic displacement jump are achieved. The (110) surface shows a somewhat more complex 
indentation behaviour; while the pertinent Schmid factor is identical to the (100) surface, the onset of plasticity rather 
follows the (111) surface. We assume that this is due to the complex stress distribution found below this surface; in 
particular, the maximum shear stress is not below the indenter along the indentation axis but rather shifted to the 
side; wc also observe more than one primary nucleation site for dislocations. Note also that it has been shown 14 ' 25 
that - in single-crystalline Cu - non-Schmid factors may be important in describing dislocation nucleation for uniaxial 
tension in [110] direction, but not in [111] or in [100] direction. 

The influence of the surface orientation on the force-displacement curves increases in the elastic regime until the 
critical stress is reached and the load drop appears. At this critical indentation depth, the orientation dependence 
is maximum, i.e., the forces exerted on the three principal surfaces for identical indentation vary most. Due to the 
stronger anisotropy of Cu, the orientation dependence is stronger for this material than for Al. Note, however, that 
also in Al, which is almost isotropic, a distinct orientation dependence is seen around the critical indentation depth; 
this effect is not due to the elastic anisotropy but rather to the processes of dislocation nucleation, see the discussion 
above. After the load drop, in both materials the force-displacement curves show almost no orientation dependence, 
apart form fluctuations, and increase linearly with depth. 

This feature is better discussed with the help of the contact pressure, that is the force divided by the projected 
contact area. Note that the measurement of this area in an atomistic simulation is not without ambiguities, see the 
Appendix. The pertinent data are displayed in Fig. 4. Initially, before the load drop, the pressure curves increase in 
agreement with the -s/d-dependence of the Hertzian theory, Eq. (3). The pressure for the (111) and the (110) surface 
increase quite similarly, in agreement with the quite similar values of the indentation moduli for these two orientations, 
cf. Table II. The pressure on the (100) surface, however, starts deviating quite early from that of the other surfaces, 
and is considerably lower; this is in agreement with the above discussion of the earlier dislocation nucleation occurring 
under this surface. The load drop shows up in a more pronounced way in this plot, since the projected areas increase 
considerably during and after the nucleation of plasticity. 

After the end of the load drop, the contact pressure reaches a rather constant value; this defines the hardness of 
the material. Thus Fig. 4 demonstrates that the surface orientation loses its influence on the hardness as soon as the 
plastic regime has been entered. The hardness of Al is measured in the simulation to be around 6 GPa, and that of Cu 
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around 12 GPa; these values correspond quite well to the simulational data obtained in Rcf. 12. Experimental data 
as obtained by nanoindentation in single-crystalline Cu (100) give only slightly smaller values, around 6 GPa, 26 while 
measurements on (ultra-fine grained) poly-crystalline specimens obtain smaller hardness values, around 2 GPa. 27 



IV. CONCLUSIONS 



We performed simulations of nanoindentation into an elastically isotropic metal, Al, and an elastically anisotropic 
metal, Cu. The three principal surfaces - the (100), (110) and (111) surface - were studied. We found: 

1. The elastic deformation is fairly well described by linear elasticity theory, as tabulated by Vlassak and Nix. 
With increasing indentation, the effect of pressure hardening is visible. 

2. Around the critical indentation depth, just before the onset of dislocation nucleation, the indentation forces 
show the strongest dependence on surface orientation. This dependence is almost equally strong in the elastic 
isotropic Al as in Cu. It is not an clastic effect but due to the fact that dislocations nucleate at different global 
stresses due to the different orientation of the slip planes. 

3. After the load drop, the contact pressure inside the material stays roughly constant and assumes the same value 
independently of the surface orientation. This demonstrates that the indentation hardness is a well defined 
concept even when the critical indentation depth has been only slightly exceeded. 
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APPENDIX A: HOW TO DETERMINE THE CONTACT AREA? 



The choice of the contact area is crucial for determining the contact pressure, but unfortunately not unique. On 
the macroscopic length scale at least two conceptionally different concepts exist: the Meyer and the Brinell contact 
areas, 3 which represent the convex contact surface and the projection onto the initial surface plane, respectively. 

For the purposes of nanoindentation, and in agreement with Hertz, the contact pressure is the substrate response 
in normal direction and therefore we have to choose the projected area A atom . On the atomistic level the obvious 
choice to measure this area is to sum up the (projected) areas of all atoms i which are in contact with the indenter, 
cf. Fig. 5: 



A atom = 7Rj 2 cosa i- ( A1 ) 



i£contact 



Here, a is an 'atom radius', and oti is the angle between the indentation direction and the vector joining the 
centre of the indenting sphere with atom i. For Cu we choose a = 1.35 A. Unfortunately, this definition leads 
to serious problems, since A atom describes a non-connected area; indeed, between the atoms there appear 'holes'. 
This becomes an important problem in particular in the elastic regime, where due to the elastic deformation, these 
'holes' increase in size with increasing deformation, leading to systematic errors: the area is underestimated, and the 
pressure is overestimated. A further systematic effect shows up due to the different areal densities of surface planes: 
For atomically rough surfaces, such as the (110), the area is underestimated in comparison to dense surfaces, such 
as the (111); this effect leads to a distortion of hardness determinations from simulation data for different surfaces. 
Fig. 6 demonstrates these effects: besides a general overestimation of the contact area, the order of the stiffness for 
the (111) and the (110) surfaces is exchanged. 

In the present paper, we therefore have chosen an elliptic contact area, which is connected and conserves the right 
order of the pressures: 



^elliptic = * ^max _ ^min^max _ y min^ ^ 
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TABLE I: Elastic properties of Al and Cu as described by the potentials by Zhou et al. and by Mishin et al. . Elastic 
constants dj, elastic anisotropy X, Eq. (5), bulk modulus B, average shear modulus G, and average Poisson ratio v. 





en (GPa) 


da (GPa) 


C44 (GPa) 


X 


B (GPa) 


G (GPa) 


V 


Al 


114 


61.5 


31.6 


1.20 


79.0 


29.3 


0.33 


Cu 


169.9 


122.6 


76.2 


3.22 


138.4 


47.8 


0.35 



Here, x and y measure the coordinates of the contact atoms projected into the initial surface plane. This set of 
coordinates describes a curved contour line, which is approximated by an ellipse. In Eq. (A2), x max — x mm and 
ymax _ denote t ne major and minor diameters of this ellipse, respectively. This definition leads initially, i.e., 
when the indenter has contact with few atoms only, to an increased noise in the determination of the area and hence 
the contact pressure, cf. Fig. 4. For larger indentations, however, it is more reliable. 
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TABLE II: Orientation dependent indentation moduli. 
M theo : after Vlassak and Nix 9 ' 10 . 

M appr : obtained according to the simplified recipe, Eq. (12). 
M MD,init. from a fit tQ Qur MD results over rf < 2 A. 

M MD ' fu11 : from a fit to our MD results over the full elastic region. 

The isotropic indentation modulus M la ° as calculated from the isotropic constants, Eqs. (9) - (11), appropriate for a polycrystal, 
is given as a reference. 





Al 


Cu 


(100) (110) (111) 


(100) (110) (111) 


M thco (Qp a ) 
M appr ( Qpa ) 


87.0 88.4 88.9 
80.5 90.3 93.4 


135.0 148.0 151.9 
81.4 179 211 


M MU - mlt (GPa) 

M MD,full (G p a) 


92 96 101 
95 97 106 


134 155 171 
147 176 192 


M' B ° (GPa) 


88.2 


145.9 





FIG. 2: Force F vs displacement d in Al (a) and Cu (b) for the three surface orientations studied. 
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d (A) d (A) 

(a) (b) 



FIG. 3: Comparison of the force F vs displacement d data as obtained by molecular dynamics (MD) with various theories and 
fits. 

theory: theoretical prediction for anisotropic media by Vlassak and Nix. 9,10 
approximate: Eq. (12). 

MD fit, initial: Fit of MD data to Hertz' law, Eq. (1), for d < 2 A. 

MD fit, full: Fit of MD data to Hertz' law, Eq. (1), for the full elastic region. 

a) Al (100). b) Cu (100). 




FIG. 5: Schematics of nanoindentation by a sphere of radius R. The indentation depth d is measured with respect to the initial 
surface plane. The relevant area A proj is the contact area of the indenter with the substrate projected into the initial surface 
plane, on is the angle between the indentation direction and the vector joining the centre of the indenting sphere with atom i. 




FIG. 6: Contact pressure vs displacement d in Cu for the three surface orientations studied, as measured with a too simplistic 
definition of the indenter contact area, Eq. (Al). 
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